C..............................<TWIEMIS.FOR>......02-JUL-1988 22:13196
C......... Copies of the command files that run this program are attached
C......... at the bottom of this file
C
 999  FORMAT(/9x,'.....................<TWIEMIS.FOR>................'
     > ,/5x,'This program uses the altitude profiles of the volume'
     > ,/5x,'emission rates from the output of  the Richards and Torr'
     > ,/5x,'interhemispheric model(the FLIP) and transforms them to'
     > ,/5x,'the line of sight from a ground station and integrates the'
     > ,/5x,'emission rates. The data must be in the form put out in'
     > ,/5x,'the FLIPRINT.COM output data files (e.g FL_GRID_N.DAT).'
     > ,/5x,'For the data to be read correctly, these files must'
     > ,/5x,'contain only one block of altitude profiles for each time'
     > ,/5x,'step. TWIEMIS.FOR was written by P Richards in January'
     > ,/5x,'1987.'/)
C
C..... SATIME = solar apparent time, RLT = SAT from FLIP
C........Xiang changed this line,enlarged the NTIM.
C........Xiang found NTIM can not be bigger than 21040.
      !PARAMETER (NTIM=1311,NALT=199,NEM=16)
      PARAMETER(NTIM=21000,NALT=199,NEM=16)
      REAL RLT(NTIM),E(NALT,NTIM,NEM),EMISS(NEM)
     > ,EMTOT(NEM),LTOSZA,SATIME,S(NALT),TIME(99),CHIN(99),ALT(NALT)
     > ,THETA(NALT),TDIFF,UTHRS
      INTEGER YEAR,DAY
      LOGICAL DEBUG
      CHARACTER*150 CHEMIS
      DATA DEBUG/.FALSE./, TDIFF/-99.0/
      DATA RE/6371.0/, ISMIN/1/
C
C....... Xiang Zhao added the following line to open input/output files
      CALL OPEN_FILE()
C....... End of Xiang's code
C....... write opening message
      WRITE(6,999)
C....... Read control data from command file
      READ(5,*) AZIM
      READ(5,*) DECL
      READ(5,*) ELEV
      READ(5,*) BEGTIM
      READ(5,*) ENDTIM
      READ(5,*) TINC
      READ(5,*) ISMAX !.. Column number of last column emission rate to output
       
C.....  Read the header then the data off the input data file which
C..... is read on unit 3. This is the ZDATAN.DAT or ZDATAS.DAT file
C..... RHEAD also prints the headers in the output files
c      OPEN(UNIT=3,STATUS='OLD',READONLY)
      CALL RHEAD(3,IDAY,GLAT,GLONG,YEAR,DAY,AZIM,ELEV,DECL,CHEMIS)

C..... Calculate number of columns of data to read
      NREM=ISMAX-1
      IF(NREM.LE.9) THEN
        IPRIN = NREM
      ELSE
        IPRIN = 9
      ENDIF
      CALL READAT(1,NTIM,NALT,NEM,RLT,ALT,E,NLTS,JMAX,NREM,TDIFF)
      WRITE(6,910) RLT(1),RLT(NLTS)
 910  FORMAT(/2X,' Emission data read in. First time=',F6.2,
     >  3X,'Last time=',F6.2/)
C
C......... make sure requested times are in the range of the data
      IF(ENDTIM.LT.BEGTIM) WRITE(6,*) ' ***** stop time < start time'
      IF(ENDTIM.LT.BEGTIM) STOP
C
      IF(ENDTIM.GT.RLT(NLTS)) WRITE(6,904) RLT(NLTS)
 904  FORMAT(/2x,'**** requested end Local Time too large, reset to
     >  max time in data =',F6.2/)
      IF(ENDTIM.GT.RLT(NLTS)) ENDTIM=RLT(NLTS)
C
      IF(BEGTIM.LT.RLT(1)) WRITE(6,909) RLT(1)
 909  FORMAT(/2X,'**** requested begin Local Time is too early, reset to
     >  min time in data=',F6.2/)
      IF(BEGTIM.LT.RLT(1)) BEGTIM=RLT(1)
C
C......... convert to radians
      GLATR=GLAT/57.269
      GLONGR=GLONG/57.269
      ELEV=ELEV/57.29578
      RAZIM=AZIM/57.29578
C
C..... The distance along the line of sight S(J)
      CALL CALOSD(BOT,TOP,DELH,ELEV,RE,JMAX,ALT,S,THETA)
C
C.......... To avoid using invalid S(J+1).
        JMAX=JMAX-1
        
C
C....... Write Headings for the files
      CALL WHEADERS(CHEMIS)
C
C.......... Loop through the requested local times: BEGTIM<= LT <=ENDTIM 
C......... at intervals TINC
      SATIME=BEGTIM-TINC
      DO 99 I=1,NTIM
      SATIME=SATIME+TINC
      IF(SATIME.GT.ENDTIM) THEN
         WRITE(6,691) BEGTIM,ENDTIM,TINC,SATIME
 691  FORMAT(X,'BEGTIM= ',F6.2,'  ENDTIM= ',F6.2,'  TINC= ',F6.2
     > ,'   SATIME= ',F6.2,/ ' SATIME > END TIME')
         STOP
      ENDIF
C..... Hour angle (HR). UTHRS determined from FLIP output
      HR=(SATIME-12.0)*15.0/57.29578
      UTHRS=TDIFF+SATIME
      CALL SOLDEC(IDAY,UTHRS,DECLR,ETRAN)
      SZA=LTOSZA(GLATR,DECLR,HR)
C....      WRITE(6,97) SATIME,UTHRS,57.3*SZA,DECLR*57.3,ETRAN,GLONG
 97   FORMAT(1X,9F9.2)
C
C
      DO 47 IS=1,NEM
      EMTOT(IS)=0.0
 47   EMISS(IS)=0.0

C
C........ Evaluate SZA, and LT along line of sight and calculate 
C........ emission rates
      DO 5 J=1,JMAX
C.......... Calculate latitude and longitude of point.
       CALL STATRA(RE,GLONGR,GLONP,DGLONG,GLATR,GLATP,ELEV,ALT(J)
     >                  ,RAZIM,S(J),THETA(J))
C..      CALL TWISUB(GLATP,GLATR,DLAT,GLONP,GLONGR,DGLONG,ALT(J)
C..     >    ,RAZIM,ELEV,J,S(J))
C..... solar zenith angle CHIN(J) and local time TIME(J) at point
C..... convert difference in long(rad) into difference in time(hr)
      HLONG=DGLONG* 3.81973
      HR=(SATIME+HLONG-12.0)*15./57.29578
      CHIN(J)=LTOSZA(GLATP,DECLR,HR)
      TIME(J)=(HR*57.29578/15.) +12.0
C
C      IF(BEGTIM.EQ.ENDTIM)  WRITE(6,771) J,S(J),GLATP,GLATR
C     > ,GLONP,GLONGR,DGLONG,CHIN(J),TIME(J)
 771    FORMAT(I3, 1P22E11.3)
C.......... ensure time is available in data
      IF(TIME(J).GT.RLT(NLTS).OR.TIME(J).LT.RLT(1)) GO TO 99
C
C............ calculate the emission rates for TIME(J) and ALT(J)
      CALL FITEM(NLTS,RLT,E,TIME(J),NTIM,NALT,NEM,J,EMISS
     >  ,ISMIN,NREM)
C...... Column emission in R
      DO 147 IS=ISMIN,NREM
      EMTOT(IS)=EMTOT(IS)+ABS(S(J)-S(J+1))*1.0E5*EMISS(IS) * 1.0E-6
 147  CONTINUE
C
      IF(ABS(BEGTIM-ENDTIM).LE.0.1*TINC) THEN
        WRITE(7,310) NINT(ALT(J)),TIME(J),CHIN(J)*57.29578
     >    ,GLONP*57.29578,GLATP*57.29578,S(J+1)
     >    ,(EMISS(IREM),IREM=1,NREM)
      ENDIF
C
 5    CONTINUE
C.......END HEIGHT LOOP
C
C
C...... Print  column emission rates 
      WRITE(9,195) SATIME+TDIFF,SATIME,SZA*57.29578,
     >   (EMTOT(IS),IS=1,NREM)

 195  FORMAT(2F8.2,F9.2,1P22E9.2)
C
       !.. calculate ratios of 6300+6364 to 7320+7330 emissions. 
       !.. Eliminated 2020-08-03
c      IF(EMTOT(1).GT.0) THEN
c        R6373=1.3*EMTOT(10)/EMTOT(1)/1.7
c        WRITE(19,195) SATIME,SZA*57.29578,R6373
c      ENDIF
C
 99   CONTINUE
 300  FORMAT(2I8,F8.2,6I8,3F10.4,1P5E10.2)
 310  FORMAT(I6,2F8.2,2F8.2,6(1P,22E9.1))
 311  FORMAT(I5,1P22E10.2)
 500  FORMAT(7F15.3)
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE CALOSD(BOT,TOP,DELH,ELEV,RE,JMAX,ALT,S,THETA)
C......... calculates the distance S(J) along the line of sight from the
C......... observer.
      LOGICAL DEBUG
      DIMENSION ALT(JMAX),S(JMAX),THETA(JMAX)
      DATA DEBUG/.FALSE./
C
      DO 3 J=1,JMAX
C.......This is the calculation of  the magnitude of vector RHO
      THETA(J)=ACOS( RE*SIN(90./57.29578+ELEV)/(RE+ALT(J)) )-ELEV
      IF (NINT(ELEV*57.29578) .EQ. 90) THEN
         S(J)=ALT(J)
      ELSE
          S(J)=(SIN(THETA(J))*(ALT(J)+RE))/SIN(90./57.29578-ELEV)
      ENDIF
      IF (DEBUG) WRITE(6,91) J,S(J),THETA(J)
 91   FORMAT(' CALOSD: ',I5,1P9E9.1)
 3    CONTINUE
      RETURN
      END
C
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      FUNCTION HANGLE ( SATIME, ELONG, EPHTR )
C......  THIS FUNCTION FINDS THE HOUR ANGLE OF THE SUN FROM APPARENT NOON.
C......  IT DOES SO BY
C......     1)  CONVERTING THE LOCAL STANDARD TIME (SATIME) TO LOCAL MEAN
C......          TIME (LMT) BY CONSIDERING THE ELONG, AND
C......     2) SUBTRACTING THE EPHEMERIS TIME ( LOCAL MEAN NOON )
C......          FROM THE LMT.
      REAL SATIME
      HSATIME=SATIME
      HEPHTR=EPHTR
      HLONG = ELONG / 15.0
      IF ( HLONG .LT. 0.0 ) HLONG = 24.0 + HLONG
      IDTOUT = IFIX(HLONG+0.50)
      ADTOUT = FLOAT(IDTOUT)
      HANGLE = 15.0 * ( HSATIME + HLONG - ADTOUT - HEPHTR )
      IF ( ABS(HANGLE) .GT. 180.0 ) HANGLE = HANGLE - SIGN(360.0,HANGLE)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      FUNCTION HRDGRS(HDMMSS)
C......     THIS FUNCTION CONVERTS HOURS ( OR DEGREES ), MINUTES AND
C......     SECONDS INTO HOURS ( OR DEGREES ).
      XMMSS = AMOD( HDMMSS, 10000.0 )
      SECS = AMOD( XMMSS, 100.0 )
      HRDGRS = ( HDMMSS - XMMSS ) / 10000.0 + ( XMMSS - SECS ) / 6000.0
     2         + SECS / 3600.0
      RETURN
      ENTRY HMSDMS(HDMMSS)
      HOURS = HDMMSS
C......     THIS FUNCTION CONVERTS HOURS ( DEGREES ) INTO HOURS ( DEGREES )
C......     MINUTES AND SECONDS
      IH = IFIX(HOURS)
      H = FLOAT(IH)
C......     THE 0.0001 ADDED BELOW IS NEEDED TO COUNTERACT ROUNDOFF ERROR
C......     IT ADDS ABOUT 0.006 SECONDS TO THE FINAL ANSWER.
      XMIN = ( HOURS - H ) * 60.0 + 0.0001
      IX = IFIX(XMIN)
      X = FLOAT(IX)
      HRDGRS = H*10000.0+X*100.0+(XMIN-X)*60.0
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      REAL FUNCTION LTOSZA(G,D,H)
C....... Function to convert LT to SZA
            LTOSZA=ACOS(SIN(G)*SIN(D)+COS(G)*COS(D)*COS(H))
      RETURN
      END
C::::::::::::::::::::::::: ZATOLT :::::::::::::::::::::::::::::::
C.......CONVERT SZA TO LT
      REAL FUNCTION ZATOLT(CHI,G,D,E)
      H=ACOS((COS(CHI)-SIN(G)*SIN(D))/(COS(G)*COS(D)))
      ZATOLT=((H-E)*57.29578/15.0)
      RETURN
      END
C::::::::::::::::::::::::::::::::: WHEADERS :::::::::::::::::::::::::::::
C..... This subroutine writes the headers in the output file
      SUBROUTINE WHEADERS(CHEMIS)
      CHARACTER*150 CHEMIS
C....... for LOS quantities file
      WRITE(7,590) CHEMIS
  590  FORMAT(//2X,'This file contains interpolated quantities at'
     > ,1X,'points along the line of sight when only one time is'
     > ,1X,'specified (FST=LST) in the TWI-RUN file. '
     > ,//,'  ALT     LT     SZA     LONG     LAT    DIST(km)',A)
      WRITE(9,495) 
 495  FORMAT(//' In the following table, columns 3 and higher are the'
     > ,/,' line of sight integrated quantities corresponding to' 
     > ,X,'columns 2 and higher from the FLIPRINT file read in.'
     > ,/,' SAT= solar apparent time (LT) and SZA= solar zenith angle'
     > ,//,1X,75('-'))

      WRITE(9,496) CHEMIS
 496  FORMAT('     UT      SAT    SZA   ',A)

C      WRITE(19,196)
 196  FORMAT(/'  Ratio of the 7320 to 6300 emissions'
     > ,/,'     SAT  SZA    6300/7320')
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE FITEM(NLTS,RLT,E,SZA,NTIM,NALT,NEM,IZ,EMISS
     >  ,ISMIN,NREM)
C......... Linear interpolation of emission rates
      REAL RLT(NTIM),EMISS(NEM),E(NALT,NTIM,NEM)
      DO 19 K=1,NLTS
      IF(SZA.GE.RLT(K).AND.SZA.LT.RLT(K+1)) GO TO 20
 19   CONTINUE
      WRITE(6,*) ' LT=',SZA, ' outside range. Emission=0.00'
      EMINT=0.0
      RETURN
C
 20   CONTINUE
      DO 30 IS=ISMIN,NREM
      A=(E(IZ,K+1,IS)-E(IZ,K,IS))/(RLT(K+1)-RLT(K))
      B=E(IZ,K+1,IS)-A*RLT(K+1)
      EMISS(IS)=A*SZA+B
 30   CONTINUE
 90   FORMAT(1P9E10.2)
      RETURN
      END 
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C......... Program to read the header from the FLIP output files
      SUBROUTINE RHEAD(ID,IDAY,GLATNR,GLONNR,YEAR,DAY,AZIM,ELEV,DECL
     >  ,CHEMIS)
      INTEGER YEAR,DAY,GLATN,GLONN,GLATS,GLONS
      CHARACTER *5 LOCAT
      CHARACTER*8 CHALT
      CHARACTER*150 CHEMIS
       
      !.. Read the headers off file and get first line of data
 991  READ(ID,902,ERR=991,END=997) YEAR,DAY,RL,F107,F107A,AP
      IF(YEAR.LT.1900.OR.YEAR.GT.4000) GOTO 991  !.. Is year OK?
 902   FORMAT(7X,I4,7X,I3,11X,F8.3,9X,F6.1,10X,F6.1,6X,F6.1)

      READ(ID,904,ERR=997,END=997) GLATN,GLONN,GLATS,GLONS
 904  FORMAT(/33X,I3,2X,I5,14X,I3,2X,I5)

C....... Now look for whether it is northern or southern hemisphere file and
C....... write it into log file
 5    READ(ID,'(6X,A5)',END=998,ERR=998) LOCAT
      IF(LOCAT.NE.'NORTH'.AND.LOCAT.NE.'SOUTH') GO TO 5
      WRITE(6,906) LOCAT 
 906  FORMAT(/3X,'******  ',A5,' emissions  **********'/)
C
      IDAY=YEAR*1000+DAY
      GLATNR=GLATN
      GLONNR=GLONN
      IF(LOCAT.EQ.'SOUTH') GLATNR=-GLATS
      IF(LOCAT.EQ.'SOUTH') GLONNR=GLONS
C
C...... solar declination
      CALL SOLDEC(IDAY,12.0,DECL,ETRAN)
      DECL=DECL*57.29578
C
C
C..... Now print headers in files
C
      WRITE(7,94) YEAR,DAY,NINT(F107),NINT(AP),RL
      WRITE(9,94) YEAR,DAY,NINT(F107),NINT(AP),RL
      WRITE(6,94) YEAR,DAY,NINT(F107),NINT(AP),RL
C      WRITE(19,94) YEAR,DAY,NINT(F107),NINT(AP),RL
 94   FORMAT(/2X,' YEAR= ',I4,',   DAY= ',I3,',   F10.7= ',I3
     > ,',  AP= ',I3,',    LSHELL= ',F6.2)
      WRITE(7,343)  NINT(ABS(GLATNR)),LOCAT,NINT(ABS(GLONNR))
     > ,NINT(AZIM),NINT(ELEV),DECL
      WRITE(6,343)  NINT(ABS(GLATNR)),LOCAT,NINT(ABS(GLONNR))
     > ,NINT(AZIM),NINT(ELEV),DECL
      WRITE(9,343)  NINT(ABS(GLATNR)),LOCAT,NINT(ABS(GLONNR))
     > ,NINT(AZIM),NINT(ELEV),DECL
C      WRITE(19,343)  NINT(ABS(GLATNR)),LOCAT,NINT(ABS(GLONNR))
C     > ,NINT(AZIM),NINT(ELEV),DECL
 343  FORMAT(/,2X,' Latitude= ',I5, X, A5 ';  Longitude= ', I5,' EAST'
     > ,//2X,' Azimuth= ',I3, ' deg,  Elevation= ',I3,' deg,'
     > ,'  solar dec= ',F7.1,' deg')

C----- look for and store the header for the data block 
      DO 420 I=1,140
        READ(ID,490,ERR=999) CHALT,CHEMIS
        IF(CHALT.EQ.'  ALT   ') THEN
           REWIND (ID)
           RETURN
        ENDIF
 420  CONTINUE
 490  FORMAT(A,A)
      RETURN
 997  WRITE(6,*) '  Error or EOF while reading data file header'
      STOP
 998  WRITE(6,*) ' ** FAILED TO FIND "NORTH" or "SOUTH" in file HEADER'
      STOP
 999  WRITE(6,*) ' ** FAILED TO FIND HEADER CONTAINING "  ALT  "'
      STOP
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE READAT(ISW,NTIM,NALT,NEM,RLT,Z,E,NLTS,MAXZ,NREM,TDIFF)
C......... Program to read in the emission rates from the files ZDATAN.DAT,
C......... ZDATAS.DAT, MDATAN.DAT and MDATAS.DAT.
      REAL UT,RLT(NTIM),Z(NALT),E(NALT,NTIM,NEM),TDIFF
      LOGICAL DEBUG/.FALSE./
C
      WRITE(6,885)
 885   FORMAT(/9X,' ********  Reading in volume emission rates ***'/)
      DO 9000 L=1,NTIM
C........ Read header to find RLT (LOCAL TIME)
      DO 9 I=1,99
      READ(3,*,ERR=9,END=9999) UT,RLT(L),CHIP
      IF(TDIFF.LT.0) TDIFF=UT-RLT(L)
      IF(L.NE.1.AND.RLT(L).LT.RLT(L-1)) THEN
 11       RLT(L)=RLT(L)+24.0
          IF(RLT(L).LT.RLT(L-1)) GO TO 11
      ENDIF
      WRITE(6,901) UT,RLT(L),NINT(CHIP)
 901  FORMAT('  UT=',F10.2,'  LT=',F6.2,'   CHI=',I5)
      GO TO 12
  9   CONTINUE
  12  CONTINUE
C
C
C........ read until find first data
      DO 19 I=1,NALT
      READ(3,*,ERR=19,END=9999) Z(1),(E(1,L,K),K=1,NREM)
      IF (DEBUG) WRITE(6,902) NINT(Z(1)),(E(1,L,K),K=1,NREM)
      GO TO 23
 19   CONTINUE
 23   CONTINUE
C
C......... Read rest of data block
      DO 39 I=2,NALT
      READ(3,*,ERR=48,END=9999) Z(I),(E(I,L,K),K=1,NREM)
      IF (DEBUG) WRITE(6,902) NINT(Z(I)),(E(I,L,K),K=1,NREM)
 39   CONTINUE
 48   CONTINUE
      IF(I.GT.MAXZ) MAXZ=I-1
      IF(ISW.EQ.0) RETURN
 9000 CONTINUE
 9999 NLTS=L-1
      IF(NLTS.GE.NTIM-1) WRITE(6,890)
 890  FORMAT(/1X,' *** The data file had too many times.'
     >  'Not all read in')
      RETURN
 902  FORMAT(I4,1P22E8.1)
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE TWISUB(GLATP,GLATR,DLAT,GLONP,GLONGR,DGLONG,ALT,
     >    RAZIM,ELEV,J,RHO)
C.........This subroutine calculates the longitudes and latitudes of the
C.........points corresponding to the altitudes in the data along the line 
C.........of sight !!!!
C.........REF : FUNDAMENTALS OF ASTRODYNAMICS  by, R.R.Bates,D.R.Mueller 
C
C
      REAL ALT,GLATP,GLATR,DLAT,GLONP,GLONGR,DGLONG,RAZIM,ELEV,
     > AE,AS,AZ,XA,XB,XC,XD,XE,XF,XG,XI,AI,AJ,AK,RHO
      LOGICAL DEBUG
      DATA DEBUG/.FALSE./
      DATA RE/6371.0/, PI/3.14159265/
C.......This part corresponds to transformation matrix 2.6_12
      XA=SIN(GLATR)*COS(GLONGR)
      XB=-SIN(GLONGR)
      XC=COS(GLATR)*COS(GLONGR)
      XD=SIN(GLATR)*SIN(GLONGR)
      XE=COS(GLONGR)
      XF=COS(GLATR)*SIN(GLONGR)
      XG=-COS(GLATR)
      XI=SIN(GLATR)
C.......Corresponds to equations 2.7_2 and 2.7_5
      AS= -RHO*COS(ELEV)*COS(RAZIM)
      AE=RHO*COS(ELEV)*SIN(RAZIM)
      AZ=RHO*SIN(ELEV)+RE
C.......Corresponds to eqn. 2.6_12 the matrice !!
      AI=XA*AS+XB*AE+XC*AZ
      AJ=XD*AS+XE*AE+XF*AZ
      AK=XG*AS+XI*AZ
C...... Component along I_K plane  i.e., latitude  of the point 
      GLATP=ATAN(AK/SQRT(AI**2+AJ**2))
C...... Component along I_J plane  i.e., longitude of the point
      GLONP=ATAN(AJ/AI)
      IF(DEBUG) WRITE(6,91) GLONP,GLONGR,GLATP,GLATR
      IF (GLONGR .EQ.PI .OR. GLONGR .GT.PI)
     > GLONP=GLONP+PI
      IF ( (RAZIM.LE.180/57.29578) .AND. ((GLONP+180./57.29578).LT.
     > 360./57.29578) )  GLONP=GLONP+180./57.29578
C....      IF(GLONP.LT.0) GLONP=PI+GLONP
      DLAT=(GLATP-GLATR)
      DGLONG=(GLONP-GLONGR)
 91   FORMAT('  TWISUB: ',1P22E11.3)
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C**STATRA**STATRA**STATRA**STATRA**STATRA**STATRA**STATRA**STATRA**STATRA
C
C
       SUBROUTINE STATRA(RE,GLONGR,GLONP,DGLONG,GLATR,GLATP,RELEV,ALTP
     >                  ,RAZIM,DISP,THETA)
C
C....       RE=6378.388
C
C...........THE GROUND DISTANCE BETWEEN THE OBSEVER'S COORDINATES AND
C...........THE GEODETIC SUB COORDINATES OF THE POINT ALONG THE LINE 
C...........OF SIGHT -- A
C
       A=THETA
C
C...........THE LATITUDE OF THE POINT ALONG THE LINE OF SIGHT
C
       B=90./57.29578-GLATR
       CC=(COS(A)*COS(B))+(SIN(A)*SIN(B)*COS(RAZIM))
       ccc=acos(cc)
       GLATP=(90./57.29578)-ccc
C
C
C...........CHANGE IN LONGITUDE
C
       DGLONG =asin(SIN(RAZIM)*SIN(A))/SIN(cCC)
C
       AZIM=RAZIM*57.263
       IF (AZIM.GE.0.0.AND.AZIM.LE.180.0) THEN
          GLONP=GLONGR+DGLONG
       ELSE
          GLONP=GLONGR-DGLONG
       ENDIF              
       GLONP=GLONGR+DGLONG
C
C...........DISTANCE ALONG THE LINE OF SIGHT
C
       DISP = SQRT(RE**2+(RE+ALTP)**2-(2*RE*(RE+ALTP)*COS(A)))
C
       RETURN
       END
C::::::::::::::::::::::::::: SOLDEC :::::::::::::::::::::::::
C...... This routine calculates the solar declination (DELTA radians)
C...... and ephemeris transit time (ETRAN in secs). ETRAN is the time
C...... the sun is overhead in Greenwich
C...... REFERENCE - page 484 "Explanatory Supplement to the Astronomical
C...... Almanac" Kenneth Seidelmann, University Science Books, 20 Edgehill 
C...... Road, Mill Valley, CA 94941, 1992
C...... IDAY is yyyyddd (eg 1988015 = feb 15), UT is the universal time
C...... in hours (0-24) 
       SUBROUTINE SOLDEC(IDAY,UT,DELTA,ETRAN)
C..... Find day of year (0-366), and year (1996), then Julian day
       INTEGER JD,DAYNUM
       REAL T,YEAR,UT,L,G,LAMBDA,EPSIL,E,GHA,DELTA,SD,ETRAN,DTOR
       DOUBLE PRECISION DJD,DUT
       DATA DTOR/57.29578/   ! degrees to radians
C...... Recover year and date and make sure UT is less than 24
       DAYNUM=MOD(IDAY,1000)
       YEAR=(IDAY-DAYNUM)/1000
       JD=INT(365.25*(YEAR-1900)+DAYNUM)+2415020      ! Julian day
       DJD=JD                                ! extra precision needed         
       DUT=UT
       T=(DJD+DUT/24.0-2451545.0)/36525.0         ! # of centuries
       L=AMOD(280.460+36000.770*T,360.0)              ! aberration
       G=AMOD(357.528+35999.050*T,360.0)              ! mean anomaly
C....... LAMBDA= ecliptic longitude. DELTA=obliquity of the ecliptic.
       LAMBDA=AMOD(L+1.915*SIN(G/DTOR)+0.020*SIN(2.0*G/DTOR),360.0)
       EPSIL=23.4393-0.01300*T
C....... Equation of time. Time difference between noon and overhead sun
       E=-1.915*SIN(G/DTOR)-0.020*SIN(2.0*G/DTOR)
     >   +2.466*SIN(2*LAMBDA/DTOR)-0.053*SIN(4*LAMBDA/DTOR)
       GHA=15*UT-180+E                    ! Greenwich hour angle
       DELTA=ASIN(SIN(EPSIL/DTOR)*SIN(LAMBDA/DTOR))  ! solar declination
       SD=0.267/(1-0.017*COS(G/DTOR))    ! 
       ETRAN=(12-E/15)*3600 ! Ephemeris transit time in secs for FLIP
       RETURN
       END